{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling Inviscid Flows with Point Sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this discussion, we'll take a look at how we can use point sources to model inviscid flows past obstacles.\n",
    "The code in this document is written in [Julia](https://www.julialang.org), but it should be straightforward to try this in any language of your choice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some Setup Code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code here is really only used to draw the pathlines.\n",
    "It might be useful to [skip](#The-Rankine-Half-Body) this section, and come back to it later if you are curious."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Install libraries if necessary\n",
    "for pkg in [\"PotentialFlow\", \"Plots\"]\n",
    "    try\n",
    "        Pkg.installed(pkg)\n",
    "    catch\n",
    "        Pkg.add(pkg)\n",
    "    end\n",
    "end\n",
    "\n",
    "using PotentialFlow \n",
    "\n",
    "# Set default plotting options\n",
    "using Plots\n",
    "default(markerstrokealpha = 0, legend = false, ratio = 1,\n",
    "        dpi = 100, size = (400, 300), grid = false)\n",
    "clibrary(:colorbrewer);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have computed the velocity field around the obstable, we will draw a couple pathlines by simulating the trajectory of passive tracer particles.\n",
    "To do this, we need to provide a couple functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "update_x!"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Given the source distribution and freestream, returns a \n",
    "function that computes the velocity of the tracer particles\n",
    "\"\"\"\n",
    "function make_state_equation(sources, U∞)\n",
    "    # Here, `x` is a vector containing the current positions\n",
    "    # of the tracer particles, and `ẋ` is vector that that\n",
    "    # their velocities will be written into.  The `t` variable\n",
    "    # stands for the current time.\n",
    "    function (ẋ, x, t)\n",
    "        fill!(ẋ, U∞)\n",
    "        induce_velocity!(ẋ, x, sources, t)\n",
    "    end\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "Given the current location and velocity of the tracers,\n",
    "compute their location after a small time step\n",
    "\"\"\"\n",
    "function update_x!(x₊, x₋, ẋ, Δt)\n",
    "    @. x₊ = x₋ + ẋ*Δt\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will then use a 4th-order Runge-Kutta scheme to simulate the trajectories of the tracers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "trace_pathlines (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function trace_pathlines(sources, U∞, start_pos, T)\n",
    "    states = [copy(start_pos) for i in 1:length(T)+1]\n",
    "    ẋs = [allocate_velocity(start_pos) for i in 1:4]\n",
    "    \n",
    "    compute_ẋ! = make_state_equation(sources, U∞)\n",
    "\n",
    "    for (i, t) in enumerate(T)\n",
    "        rk4!(states[i+1], states[i], t, step(T), compute_ẋ!, update_x!, ẋs)\n",
    "    end\n",
    "    \n",
    "    transpose(hcat(states...))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "component"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Compute the component of `w` in the `n̂` direction\n",
    "\"\"\"\n",
    "component(n̂, w) = real(conj(n̂)*w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Rankine Half Body"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\newcommand{\\compose}[2]{#1 \\circ #2}\n",
    "\\newcommand{\\ddp}[2]{\\frac{\\partial #1}{\\partial #2}}\n",
    "\\newcommand{\\dpt}[1]{\\ddp{#1}{t}}\n",
    "\\newcommand{\\dd}[2]{\\frac{\\mathrm{d}#1}{\\mathrm{d}#2}}\n",
    "\\newcommand{\\ddt}[1]{\\dd{#1}{t}}\n",
    "\\renewcommand{\\vec}[1]{\\boldsymbol{#1}}\n",
    "\\newcommand{\\cross}{\\times}\n",
    "\\newcommand{\\abs}[1]{\\left|#1\\right|}\n",
    "\\renewcommand{\\d}{\\,\\mathrm{d}}\n",
    "\\renewcommand{\\div}{\\nabla \\cdot}\n",
    "\\renewcommand{\\curl}{\\nabla \\cross}\n",
    "\\newcommand{\\grad}{\\nabla}\n",
    "\\newcommand{\\lap}{\\nabla^2}\n",
    "\\newcommand{\\uvec}[1]{\\vec{#1}}\n",
    "\\newcommand{\\utangent}{\\uvec{\\tau}}\n",
    "\\newcommand{\\unormal}{\\uvec{n}}\n",
    "\\newcommand{\\im}{\\mathrm{i}}\n",
    "\\newcommand{\\eu}{\\mathrm{e}}\n",
    "\\newcommand{\\pint}{\\int}\n",
    "\\newcommand{\\inv}[1]{{#1}^{-1}}\n",
    "\\newcommand{\\conj}[1]{#1^\\star}\n",
    "\\newcommand{\\Res}[2]{\\mathrm{Res}\\left(#1,#2\\right)}\n",
    "\\newcommand{\\real}[1]{\\mathrm{Re}\\left\\{#1\\right\\}}\n",
    "\\newcommand{\\imag}[1]{\\mathrm{Im}\\left\\{#1\\right\\}}\n",
    "\\newcommand{\\normal}[1]{#1^\\perp}\n",
    "\\newcommand{\\tangent}[1]{#1^\\parallel}\n",
    "\\newcommand{\\residue}[2]{\\mathrm{Res}\\left(#1, #2\\right)}\n",
    "$$\n",
    "Let's start with a simple freestream in the $x$-direction with velocity $U_\\infty$, and a point source at the origin.\n",
    "What should the strength of this point source be in order to make $(-1, 0)$ a stagnation point?\n",
    "\n",
    "Recall that the velocity field induced by a point source at $\\vec{x}_s$ with strength $S$ is given by:\n",
    "\n",
    "$$\n",
    "    \\vec{u}_S(\\vec{x})\n",
    "  = \\frac{S}{2\\pi}\n",
    "    \\frac{\\vec{x} - \\vec{x}_s}\n",
    "         {\\abs{\\vec{x} - \\vec{x}_s}}.\n",
    "$$\n",
    "\n",
    "We want the velocity induced by the freestream and the point vortex to be zero at $(-1, 0)$, i.e.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}U_\\infty \\\\ 0\\end{bmatrix} + \\vec{u}_s\\left(\\begin{bmatrix}-1 \\\\ 0\\end{bmatrix}\\right) = 0\n",
    "$$\n",
    "\n",
    "Since $\\vec{x}_s = \\begin{bmatrix} 0 & 0 \\end{bmatrix}^T$, we have\n",
    "\n",
    "$$\n",
    "U_\\infty - \\frac{S}{2\\pi} = 0,\n",
    "$$\n",
    "\n",
    "which give us $S = 2\\pi U_\\infty$.\n",
    "Let's see what that looks like.\n",
    "\n",
    "We first set the freestream velocity and create a single point source:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Source.Point(0.0 + 0.0im, 6.283185307179586)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U∞ = 1.0\n",
    "source = Source.Point(0, 2π*U∞)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will seed 30 tracer particles along the $x=-3$ line, spanning $y \\in [-3, 3]$, and track their trajectory over 20 seconds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tracer_start = -3.0 .+ linspace(-3,3,30)*im;\n",
    "T = 0:1e-2:20\n",
    "\n",
    "freestream = trace_pathlines((), U∞, tracer_start, T)\n",
    "freestream_with_source = trace_pathlines(source, U∞, tracer_start, T);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt = plot(layout = 2, size = (600, 300))\n",
    "plot!(plt[1], freestream, color = colorant\"lightgray\", xlabel = \"x\", ylabel = \"y\",\n",
    "      xlim = (-3, 3), ylim = (-3, 3), title = \"Freestream\")\n",
    "plot!(plt[2], freestream_with_source, color = colorant\"lightgray\", xlabel = \"x\", ylabel = \"y\",\n",
    "      xlim = (-3, 3), ylim = (-3, 3), title = \"Freestream + Source\")\n",
    "plot!(plt[2], [source], markersize = 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What would have happened if the point source was placed off the $y$-axis?\n",
    "We will then have\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "U_\\infty - \\frac{S}{2\\pi} & = 0 \\\\\n",
    "\\frac{S}{2\\pi} \\frac{y}{1 + y^2} & = 0.\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "This means that we cannot zero out the horizontal velocity without imparting a vertical velocity.\n",
    "So instead of a stagnation point, we have a point that stops flow in only one direction.\n",
    "That's not really useful on its own, but in the next section we will see how we can use whole bunch of point sources to stop flow from going through a whole bunch places."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A whole bunch of point sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we want to model a solid surface in an inviscid fluid.\n",
    "The only role of that surface is to prevent fluid from going through it.\n",
    "We can accomplish the same thing by replacing the surface with a **source sheet**.\n",
    "We haven't covered source sheets yet, but you can think of it as a source where the strength is distributed along a curve instead of being concentrated at a single point.\n",
    "Since we will be discretizing the sheet into small segments, and represent those segments as point sources anyways, you can temporarily forget what I said about the sheets and just think of it as a lot of point sources."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by considering flow past a equilateral triangle.\n",
    "We will discretize this triangle into small segments, each with length $\\Delta s$.\n",
    "This will look something like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Δs = 0.1\n",
    "s = -(1 - Δs):Δs:(1 - Δs)\n",
    "\n",
    "bottom = Source.Point.(s, 1.0);\n",
    "right  = Source.Point.((s + 1).*exp(im*2π/3) + 1, 1.0);\n",
    "left   = Source.Point.((s .+ 1)*exp(-im*2π/3) + im*√3, 1.0);\n",
    "\n",
    "sources = vcat(bottom, right, left)\n",
    "targets = Elements.position.(sources);\n",
    "\n",
    "n̂ = vcat(\n",
    "    fill(-1.0im, length(bottom)),\n",
    "    fill(exp(im*π/6), length(right)),\n",
    "    fill(exp(im*5π/6), length(left)));\n",
    "\n",
    "Δs = fill(Δs, length(n̂));\n",
    "\n",
    "plot(sources, ratio = 1)\n",
    "quiver!(real.(targets), imag.(targets), quiver = 0.15.*(real.(n̂), imag.(n̂)), color = :black)\n",
    "plot!(sources, markersize = 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where the red cross marks correspond to the center of each segment (which we represent with a single point source), and the arrows correspond to the normal vectors along the surface.\n",
    "\n",
    "So the velocity at any point in the fluid is equal to the freestream plus the velocity induced by each of the $N$ point sources:\n",
    "\n",
    "$$\n",
    "\\vec{u}(\\vec{x}) = \\vec{U}_\\infty + \\sum_{n = 1}^N \\vec{u}_n(\\vec{x})\n",
    "$$\n",
    "\n",
    "If we want enforce the no-flow-through condition at the point sources, then we must have\n",
    "\n",
    "$$\n",
    "\\sum_{n = 1}^N \\uvec{n}_i \\cdot \\vec{u}_n(\\vec{x}_i) = - \\uvec{n}_i \\cdot \\vec{U}_\\infty \\qquad \\text{for } i = 1, 2, \\cdots, N.\n",
    "$$\n",
    "\n",
    "As long as the we are not evaluating the velocity that a segment induces on itself, we will treat the segment as a point source:\n",
    "\n",
    "$$\n",
    "    \\vec{u}_n(\\vec{x})\n",
    "  = S_n \\cdot \\underbrace{\\frac{1}{2\\pi}\n",
    "    \\frac{\\vec{x} - \\vec{x}_n}\n",
    "         {\\abs{\\vec{x} - \\vec{x}_n}}}_{\\uvec{u}_n(\\vec{x})} \\qquad \\text{when } \\vec{x} \\neq \\vec{x}_n\n",
    "$$\n",
    "\n",
    "But how **would** we evaluate $\\vec{u}_n(\\vec{x}_n)$?  We probably won't go over it in class, so you'll just have to believe ([or not](#Or-Not)) that it somehow is:\n",
    "\n",
    "$$\n",
    "\\vec{u}_n(\\vec{x}_n) = \\frac{S_n}{\\Delta s}\n",
    "$$\n",
    "\n",
    "<a id=\"back\"></a>\n",
    "\n",
    "So the no-flow-through condition above can be written in matrix form:\n",
    "\n",
    "$$\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1/\\Delta s & \\uvec{n}_1 \\cdot \\frac{1}{2\\pi}\\frac{\\vec{x}_1 - \\vec{x}_2}{\\abs{\\vec{x}_1 - \\vec{x}_2}^2}\n",
    "& \\cdots & \\uvec{n}_1 \\cdot \\frac{1}{2\\pi}\\frac{\\vec{x}_1 - \\vec{x}_N}{\\abs{\\vec{x}_1 - \\vec{x}_N}^2}\n",
    "\\\\\n",
    "\\uvec{n}_2 \\cdot \\frac{1}{2\\pi}\\frac{\\vec{x}_2 - \\vec{x}_1}{\\abs{\\vec{x}_2 - \\vec{x}_1}^2} & 1/\\Delta s\n",
    "& & \\vdots\n",
    "\\\\\n",
    "\\vdots & & \\ddots & \\vdots \\\\\n",
    "\\uvec{n}_N \\cdot \\frac{1}{2\\pi}\\frac{\\vec{x}_N - \\vec{x}_1}{\\abs{\\vec{x}_N - \\vec{x}_1}^2} &\n",
    "\\cdots & \\cdots &\n",
    "1 / \\Delta s\n",
    "\\end{bmatrix}}_{A}\n",
    "\\begin{bmatrix}\n",
    "S_1 \\\\\n",
    "S_2 \\\\\n",
    "\\vdots \\\\\n",
    "\\vdots \\\\\n",
    "S_n\n",
    "\\end{bmatrix} =\n",
    "\\underbrace{ -\n",
    "\\begin{bmatrix}\n",
    "\\uvec{n}_1 \\cdot \\vec{U}_\\infty \\\\\n",
    "\\uvec{n}_2 \\cdot \\vec{U}_\\infty \\\\\n",
    "\\vdots \\\\\n",
    "\\vdots \\\\\n",
    "\\uvec{n}_N \\cdot \\vec{U}_\\infty \\\\\n",
    "\\end{bmatrix}}_{b},\n",
    "$$\n",
    "\n",
    "which we can solve to find the strengths of the point sources.\n",
    "In code, this will look something like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b = -component.(n̂, U∞);\n",
    "\n",
    "A = [component(n, induce_velocity(target, source, 0)) for (n, target) in zip(n̂, targets), source in sources]\n",
    "A .+= 0.5*diagm(1./Δs);\n",
    "\n",
    "S = A \\ b;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sources = Source.Point.(targets, S)\n",
    "\n",
    "plot(sources, color = :RdBu_r, colorbar = :right, markersize = 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the figure above, red x's correspond to sources while blue x's correspond to sinks.\n",
    "Now, we will seed a line of tracer particles along $x = -3$ and advect them in time to trace out pathlines:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tracers = trace_pathlines(sources, U∞, tracer_start, T)\n",
    "\n",
    "plot(tracers, color = colorant\"lightgray\", xlabel = \"x\", ylabel = \"y\",\n",
    "     xlim = (-3, 3), ylim = (-3, 3), size = (500, 300))\n",
    "plot!(sources, markersize = 5, color = :RdBu_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Food for thought**\n",
    ">\n",
    "> - Did the sources *have* to be on the surface of the body?\n",
    "> - Can we achieve the same result with $N$ point sources placed *inside* the body?\n",
    "> - Why not place the sources outside the body?\n",
    "\n",
    "Before moving on, we'll wrap all of this into a function so we don't have to keep repeating ourselves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "simulate_flow (generic function with 1 method)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function simulate_flow(sources, Δs, n̂, U∞; tracer_start = -3.0 .+ linspace(-3,3,30)*im, T = 0:1e-2:20)\n",
    "    b = -component.(n̂, U∞);\n",
    "\n",
    "    A = [component(n, induce_velocity(target, source, 0)) for (n, target) in zip(n̂, targets), source in sources]\n",
    "    A .+= 0.5*diagm(1./Δs);\n",
    "\n",
    "    S = A \\ b;\n",
    "    \n",
    "    sources = Source.Point.(targets, S)\n",
    "    \n",
    "    tracers = trace_pathlines(sources, U∞, tracer_start, T)\n",
    "\n",
    "    plt = plot(tracers, color = colorant\"lightgray\", xlabel = \"x\", ylabel = \"y\",\n",
    "               xlim = (-3, 3), ylim = (-3, 3), size = (500, 300))\n",
    "    plot!(plt, sources, markersize = 2, color = :RdBu_r)\n",
    "    \n",
    "    plt, sources, tracers\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow Past a Cylinder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we first construct the surface of the body with point sources:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "θ = linspace(0, 2π, 101)[1:end-1]\n",
    "Δs = fill(step(θ), length(θ))\n",
    "sources = Source.Point.(exp.(im.*θ), 1.0);\n",
    "targets = Elements.position.(sources)\n",
    "n̂ = Elements.position.(sources);\n",
    "\n",
    "plot(sources, ratio = 1)\n",
    "quiver!(real.(targets), imag.(targets), quiver = 0.15.*(real.(n̂), imag.(n̂)), color = :black)\n",
    "plot!(sources, markersize = 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now with a freestream from left to right, we have;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt, actual_sources, tracers = simulate_flow(sources, Δs, n̂, 1.0)\n",
    "plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now with a downward freestream, this becomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt, actual_sources, tracers = simulate_flow(sources, Δs, n̂, -im, tracer_start = linspace(-3,3,40) .+ 3im)\n",
    "plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow Past a Happy Face"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "θ = linspace(π/6, 5π/6, 100)\n",
    "Δs = step(θ)\n",
    "targets = exp.(-im*θ)\n",
    "n̂ = targets;\n",
    "targets = vcat(targets, flipdim((1 - Δs)*targets, 1))\n",
    "n̂ = vcat(n̂, -flipdim(n̂, 1))\n",
    "n̂[1] *= im\n",
    "n̂[50] *= -im\n",
    "\n",
    "Δs = abs.(gradient(targets))\n",
    "\n",
    "θᵢ = linspace(0, 2π, 30)[1:end-1]\n",
    "left_eye = 0.1exp.(im*θᵢ) - 0.5 + 0.2im\n",
    "right_eye = 0.1exp.(im*θᵢ) + 0.5 + 0.2im;\n",
    "\n",
    "append!(targets, left_eye)\n",
    "append!(targets, right_eye)\n",
    "append!(Δs, fill(0.5π/length(θᵢ), 2length(θᵢ)))\n",
    "append!(n̂, exp.(im*θᵢ))\n",
    "append!(n̂, exp.(im*θᵢ))\n",
    "\n",
    "sources = Source.Point.(targets, 1.0);\n",
    "plot(sources, ratio = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With $U_\\infty$ going downward:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt, actual_sources, tracers = simulate_flow(sources, Δs, n̂, -1.0im, tracer_start = linspace(-3,3,41) .+ 3im)\n",
    "plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With $U_\\infty$ going from left to right:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt, actual_sources, tracers = simulate_flow(sources, Δs, n̂, 1.0)\n",
    "plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Appendix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Or Not"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll have to discuss source sheets for a little bit here to find out the velocity that the sheet induces on itself.\n",
    "To keep it simple, we will use a straight sheet segment with length $\\Delta s$ centered around the origin.\n",
    "The equation describing this sheet segment is given by:\n",
    "\n",
    "$$\n",
    "\\vec{u}_s(\\vec{x}) = \\frac{1}{2\\pi}\n",
    "\\int_{-\\Delta s/2}^{\\Delta s/2}\n",
    "\\sigma(l) \\frac{\\vec{x} - \\vec{x}_s(l)}\n",
    "{\\abs{\\vec{x} - \\vec{x}_s(l)}}\n",
    "\\d{l},\n",
    "$$\n",
    "\n",
    "where $\\sigma$ is the local strength of the sheet.\n",
    "This is not all that different from the point sources, except when you approach the sheet.\n",
    "\n",
    "We first assume that $\\Delta s$ is small enough that $\\sigma$ is essentially constant in the interval.\n",
    "Then, if we parametrize the target point as\n",
    "\n",
    "$$\n",
    "\\vec{x} = r\\cos\\theta\\,\\uvec{i} + r\\sin\\theta\\,\\uvec{j},\n",
    "$$\n",
    "\n",
    "we obtain\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\vec{u}_s(\\vec{x})\n",
    "& = \\frac{\\sigma(0)}{2\\pi}\n",
    "\\left[ \\uvec{i}\n",
    "\\int_{-\\Delta s/2}^{\\Delta s/2}\n",
    "\\frac{(r\\cos\\theta - l)}\n",
    "{\\abs{\\vec{x} - \\vec{x}_s(l)}}\n",
    "\\d{l} + \\uvec{j}\n",
    "\\int_{-\\Delta s/2}^{\\Delta s/2}\n",
    "\\frac{r\\cos\\theta}\n",
    "{\\abs{\\vec{x} - \\vec{x}_s(l)}}\n",
    "\\d{l}\n",
    "\\right] \\\\\n",
    "& =-\\frac{\\sigma(0)}{2\\pi}\n",
    "\\left( \\frac{\\uvec{i}}{2}\n",
    "\\log\\left[ (r\\cos\\theta - l)^2 + (r\\sin\\theta)^2\n",
    "\\right]_{-\\Delta s/2}^{\\Delta s/2} + \\left. \\uvec{j}\n",
    "\\tan^{-1 }\\left(\\frac{r\\cos\\theta - l}\n",
    "{r\\sin\\theta}\n",
    "\\right)\n",
    "\\right|_{-\\Delta s/2}^{\\Delta s/2}\n",
    "\\right)\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "As $\\vec{x}$ approaches the origin, $r$ goes to zero\n",
    "\n",
    "$$\n",
    "\\lim_{r \\rightarrow 0}\n",
    "\\vec{u}_s(\\vec{x}) = \\mathrm{sign}(\\theta)\\frac{\\sigma(0)}{2}.\n",
    "$$\n",
    "\n",
    "Since we are using point sources to model a segment of a source sheet, its strength $S$ is distributed along the segment, i.e.\n",
    "\n",
    "$$\n",
    "S = \\int_{-\\Delta s/2}^{\\Delta s/2} \\sigma(l) \\d{l}.\n",
    "$$\n",
    "\n",
    "For a short sheet segment, we assume that $\\sigma$ is essentially constant along the segment, so we have $\\sigma \\approx \\frac{S}{\\Delta s}$, which finally gives us\n",
    "\n",
    "\\begin{equation}\n",
    "\\lim_{r \\rightarrow 0} \\uvec{u}_s(\\vec{x}) = \\frac{S}{2\\Delta s}\n",
    "\\end{equation}\n",
    "\n",
    "You can go [back](#back) now."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0-pre.beta",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
